Transfer-matrix renormalization group study of the spin ladders 
with cyclic four-spin interactions 
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The temperature dependence of the specific heat and spin susceptibility of the spin ladders with 
cyclic four-spin interactions in the rung-singlet phase is explored by making use of the transfer- 
matrix renormalization group method. The values of spin gap are extracted from the specific heat 
and susceptibility, respectively. It is found that for different relative strength between interchain 
and intrachain interactions, the spin gap is approximately linear with the cyclic four-spin interaction 
in the region far away from the critical point. Furthermore, we show that the dispersion for the 
one-triplet magnon branch can be obtained by numerically fitting on the partition function. 



PACS numbers: 75.10.Jm, 75.40.Cx, 75.40.Mg 
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I. INTRODUCTION 



In the past few years, cyclic spin exchanges in quantum 
antifcrromagnctic magnets have attracted considerable 
attention. From various experiments, it has been realized 
that four-spin cyclic exchanges play an important role in 
understanding the spin dynamics of insulating cupratc 
materials, such as the two-leg spin ladders CU2O3 in the 
compounds SrC^C^? 1 - (Sr, Ca, La)i4Cu2404i f 2 -^ and 
the CuC>2 plane in La2CuC>4.— & Similar mechanism of 
multiple spin exchanges in 3 He (Ref0) and Wigner crys- 
tal^ which was found to be large, has also been proposed. 

Spin ladders can be viewed as intermediates between 
one-dimensional spin chains and two-dimensional spin 
lattices. They arc particular of interest in the under- 
standing of the high-temperature superconductivity in 
two dimensions and the Haldanc's conjcctur o 9 ' 10 in spin 
chains. The two-leg spin ladder is the simplest sys- 
tem on which cyclic exchanges can be materialized. In 
the absence of cyclic exchanges, with nonzero interchain 
couplings, the 5=1/2 spin ladder lies in a Haldane 
phase, which has a short-range resonating-valence-bond 
(RVB) ground state with a spin gap, and superconduc- 
tivity emerges upon doping. 11 ' 12 Contrary to the case 
in a 5 = 1/2 chain, the spinous in the ladder are con- 
fined back to the well-defined quasiparticlcs (massive 
magnons) in the neighborhood of q — n, where a spec- 
tral gap exists^ It has been shown that four-spin cyclic 
exchanges can make the spin gap decreases rapidly^ Fur- 
thermore, tuning of the four-spin interactions can drive 
the system into different phases! 14 ' 15 At the nearest crit- 
ical point, the system becomes gapless to cross from the 
RVB rung-singlet phase to a staggered dimer phase. By 
a field theory using the Majorana Fermion representa- 
tion, Nersesyan and Tsveliki^ referred to the critical phe- 
nomenon as the level k = 2 SU(2) Wess-Zumino-Novikov- 
Witten model or SU(2) c = | conformal field theory. 
In other words, the phase transition is of the second- 
order type and belongs to the universality class of the 
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FIG. 1: Schematic structure of a two-leg ladder with cyclic 
(four-spin) interaction. The round arrows stand for the ac- 
tion of the four-spin terms as given in Eq. (fl]), which is not 
completely identical to the usual cyclic permutations. 



Takhtajan-Babujian critical point^ This statement has 
been confirmed by numerical simulations! 18 ' 19 

The Hamiltonian describing an isotropic 5=1/2 spin 
ladder with additional cyclic (four-spin) interaction is 
represented by 



H 



{JrungSl,i ' S2,i + ^lcg [Si,, ■ Sl^+l + S24 ■ S2,i+1 



2J C yc [(Sl,i ' Si^+i) (S2,i • S2,i:+l) 
(Si,, • S2,i) (Sl,i+1 • S2,j+l) 
(Sl,i • S2,i+l) (S2,i ■ Sl,i+l)]} , 



(1) 



where the indices 1 and 2 distinguish upper and lower 
legs, and i labels rungs; J rung and Ji eg are the interchain 
and intrachain couplings, respectively; J cyc denotes the 
strength of the cyclic (four-spin) exchange. The model is 
schematically shown in Fig. [TJ In the rest of the paper, 
for convenience, we call Eq. {!]) as type-I Hamiltonian. 

The above Hamiltonian distinguishes from another one 
where the ring exchange term is directly represented by 
a cyclic permutation operator Pijki in the plaqucttc: 



K P, 



ijkl 



P~ L - - 
ijkl a 



2 



We call this form as type-II Hamiltonian. The cyclic per- 
mutation operator can be expressed by spin-1 operators 
including bilinear frustrating and biquadratic terms^^ 
So type-II Hamiltonian can be re-written as 

H = ^2 {(^rung + 2K) Si,j • S 2 ,i 

i 

+ (Jlcg + K) (Si^ • Si ;! ; + i + S2,i • S2,i+l) 

+ K (Si,i ■ S2 : i+i + S 2j i ■ Si,i+i) 

+ 4^[(s M -s M+ i)(s 2ii -s 2 , i+ i) 

+ (Sl,i • S2,i) (Si ;2 :+i • S 2 ,i+l) 

— (Si,, • S2,i+l) (S2,i ' Si^+i)]} . (2) 

There have been extensive theoretical studies on the 
type-II Hamiltonian. By the density matrix renormaliza- 
tion group (DMRG)2i and exact diagonalization method, 
the energy spectra about one-magnon branch and two- 
magnon continuum in the rung-singlet phase have been 
obtainedi^^ The phase diagram for J rung = J\ cg was 
studied in detailed 1 ^ and a global picture for general cases 
has been proposed^ The phase transitions, mainly con- 
cerning the case from the rung-singlet to the dimerized 
phase, together with the critical exponent were stud- 
ied by field theory, perturbative and series expansion 
method.i^^^ 

On the other hand, the type- 1 Hamiltonian catches less 
attention. And it is widely believed that the physics 
should be the same for the type-I and type-II Hamil- 
tonians, although there is a lack of detailed compari- 
sions. Theoretically, the ring exchange interaction natu- 
rally arises as an effective low-energy model for the multi- 
band Hubbard model near the insulating filling and in 
the strong coupling limit. It has been shown that the 
four-spin terms are most significant onesi 2 ^ In order to 
highlight the effect of this four-spin term, in this paper, 
we employ the type-I Hamiltonian for our numerical sim- 
ulations. 

As mentioned above, the zero-temperature behaviors 
of the spin ladders with cyclic exchange interactions have 
been analyzed in detail, especially for the typc-II Hamil- 
tonian. A study on the thermodynamic properties of 
the systems has not been carried out extensively, ex- 
cept for the study by a high-temperature series expan- 
sion method^ These quantities can be directly com- 
pared with the experimental observations and therefore 
are of particular use. Detailed and reliable numerical re- 
sults on these quantities may offer an expedient way to 
determine the model parameters for the compounds we 
are interested in. In this paper, we show that by using the 
transfer-matrix renormalization group (TMRG) method, 
we can not only obtain the temperature dependence of 
the specific heat and susceptibility of the spin ladders 
down to the low-temperature regime and with high ac- 
curacy, but also extract extra information from the data 
for the magnitude of the gap and the one-triplet magnon 
excitations, which can be compared with the experiments 
directly. 



As proposed by experiments, for the insulator cupratc 
ladder materials, h 2 < 3 < 4 i 25 both J mng and Ji eg are positive 
and close in magnitude. The values of K/J lung range in 
a region 0.025 < K/J Iung < 0.075, which roughly corre- 
sponds to the parameter region 0.05 < Jcyc/>7rung < 0.15 
for the type-I Hamiltonian and the system is still in the 
rung-singlet phase. For this reason, in our numerical sim- 
ulations, we confine the parameter x cyc £ [0,0.1], and x 
is taken three different values: 0.5, 1.0 and 1.2, where 

X . J\cg/ ^rungj (3) 
*£cyc ■ — J eye/ Jrung- (4) 

We take x = 0.5 in order to make a comparison with the 
typical rung-singlet phase. 

The numerical algorithm we employ to study the ther- 
modynamical properties of the present model is the 
TMRG method ) 26 i 27 ' 28 which is a powerful numerical 
tool for studying the thermodynamic properties of one- 
dimensional quantum systems. One is able to evalu- 
ate nearly all thermodynamic quantities by the max- 
imum eigenvalue and the corresponding left and right 
eigenvectors of the transfer matrix. We will show later 
that the magnitude of the gap and the low-lying excita- 
tions can also be extracted from the obtained data. In 
our numerical calculations, t = 0.05, the error caused 
by the Trotter-Suzuki decomposition is less than 10~ 3 . 
During the TMRG iterations, 80 — 100 states are re- 
tained and the truncation error is less than 10 -4 down 
to k B T ~ 0.01J rung . 



II. SPECIFIC HEAT AND SPIN 
SUSCEPTIBILITY 

Figure [2] shows the TMRG results on the temperature 
dependence of the specific heat coefficient C/T for x = 
0.5, 1.0, 1.2 with x cyc = 0, 0.02, 0.04, 0.06, 0.08 and 
0.1. Besides the exponential decrease as T — ► 0, which 
indicates the existence of a gap, the general behavior of 
the curves for a given x is as follows: As x cyc increased, 
the slope of (C/T) max shifts to lower temperatures. The 
reason is obvious: that with the increase of the four-spin 
cyclic exchange, the spin gap and the whole dispersion 
decreasei 2 ^ Another interesting feature is that at x = 
1.0 and 1.2, a shoulder gets more and more distinct at 
low temperatures with the increase of x cyc . We attribute 
this shoulder to the tendency of deconfmcmcnt of bound 
spinous from the magnons, as the crossover occurs from 
the strong coupling regime (Ji eg /"/rung = iCl) to the 
weak coupling one (</i og / Jrung — x 1) with the extreme 
limiting case of a finite C/T value at T — 0. 

The corresponding results for spin susceptibility x(T) 
are presented in Fig. [31 The exponential decrease of x(T) 
as T — > 0, and the shift of the peak to lower tempera- 
tures with the increase of consistent with the 
above behaviors of C /T. Our numerical result on the 
susceptibility well agrees with that of Ref. [24| obtained 
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FIG. 2: (Color online) The specific heat coefficient C/T for 
type-I Hamiltonian for three values of x with x cyc = 0, 0.02, 



0.04, 0.06, 0.08 
of arrows. 



and 0.1, in ascending order along the direction 



FIG. 4: 1/T dependence of In (CT 3/2 ) and In (xT 1/2 ^ (in 

the unite of J rung ) for x = 0.5 with x cyc = 0, 0.02, 0.04, 0.06, 
0.08, and 0.1 in ascending order along the direction of arrows. 



by high-temperature series expansion and exact diago- 
nalization. 

The Xcyc dependence of the gap is a central quantity, 
since basically it determines the relevancy of the cyclic 
exchange term. In the vicinity of the critical point x^. yc , 
generally, 
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FIG. 3: Temperature dependence of the susceptibility \ f° r 
type-I Hamiltonian for three values of x with x cyc = 0, 0.02, 
0.04, 0.06, 0.08, and 0.1 in ascending order along the direction 
of arrows. 



For the type-II Hamiltonian in both strong and weak cou- 
pling limits, by the perturbative cluster expansion and 
Pade approximation in the strong coupling regime, to- 
gether with the bosonization technique in the weak cou- 
pling limit, Miiller et al£2- showed that the critical expo- 
nent r\ for the gap is approximately equal to 1; i.e., the 
gap is linearly dependent on x cyc near the critical point. 
For the J\ cg = J rU ng case, this result was also obtained 
by Hijii et alr^ Here, we will see that for type-I Hamilto- 
nian, in the region considerably far away from the critical 
point, x C yc G (0,0.1), the linear behavior of the gap vs. 



L cyc 



still holds in almost all the parameter range. 



For the spin-i ladder in the rung-singlet phase, sim- 
ilar to the case in a S = 1 chain, we can assume that 
the low-lying excitation spectrum e(k) for k near 7r is 
approximately given by 



e{k) = ^Jv 2 {k-n) 2 +A 2 



<D(\k-n\ 



(5) 



where v is the spin wave velocity and A is the spectrum 
gap. The low temperature behaviors of the specific heat 
and susceptibility per spin can be expressed, respectively, 
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C(T) 



3A / A 



2irT 



3/2 



-A/T 



-A/T 



(6) 
(7) 



for T < A. In this situation, In (CT 3 / 2 ) and In (xT 1 / 2 ) 
should be linear to 1/T with a slope — A. Figure [4] shows 
the 1/T dependence of In (CT 3 / 2 ) and In (xT 1 / 2 ) (in 
units of J rU ng) f° r % = 0.5 with a;, 
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0, 0.02, 0.04, 

0.06, 0.08, and 0.1. We observe a well-behaved linearity 
at 1/T > 10. For x = 1.0 and 1.2, we have the simi- 
lar results. The values of the gap at various x and x cyc 
can thus be extracted by fitting the data in the interval 
(3 = 1/T e (10, 20). The results for x = 0.5 and x = 1.0 
are presented in Fig. [5] The gaps obtained from the spe- 
cific heat and susceptibility are well coincident with each 
other. 

To both verify the low-temperature scaling behavior of 
the system and check the validity of Eq. (J7|) , in Fig. [6j the 
rescalcd susceptibility x/Xmax versus T/A for x = 0.5, 
1.0 at various x cyc is presented. For a fixed value of 
x, the curves converge to a single one at low tempera- 
tures, demonstrating the dominance of A in determining 
the susceptibility at low temperatures. Furthermore, it 
also suggests that the values of the gap we obtained by 
this approach are quantitatively reliable. This situation 
is quite similar to the case of the S = 1 antifcrromag- 
netic Hcisenberg chain with single-ion anisotropy^ Here 
Xmax serves as a common factor which incorporates the 
unknown effect on the prcfactor in Eq. (|7|). 

From Fig. [3J we observe that the gap is approxi- 
mately linearly dependent on x cyc for fixed x in the pa- 
rameter space we studied. We also notice the differ- 
ent descending rates of A with x cyc for various x. At 
^cyc = 0, A(x = 0.5) > A(x = 1.0) with the difference 
~ 0.15. This indicates that the gap decreases as the in- 
trachain coupling J\ cg increases, which is consistent with 
the perturbative picture from the dimer limit of uncou- 
pled rungs (strong coupling regime). On the other hand, 
at x cyc = 0.1, the difference of the gaps becomes very 
small. We expect for higher value of a; cyc a reverse of the 
gaps with A(x = 0.5) < A(x = 1.0) should occur. 

For x = 1.0 and 1.2 (not shown in Fig. [5]), the values 
of the gap obtained by the present method well coincide 
with the numerical results from the density-matrix renor- 
malization group (DMRG) study for the typc-II Hamilto- 
nian for the corresponding values of the cyclic exchange 
couplingjiiii For instance, when x = 1.0 and x cyc = 0, 
the gap obtained here by TMRG is ~ 0.51, as compared 
with the DMRG result of 0.504.±i 

To come to a close, we estimates the errors in the pro- 
cess of gap fitting. On the one hand, as mentioned in the 
previous section, combining the errors from the Trotter- 
Suzuki decomposition with the TMRG truncation, we 
estimate that the error caused by the TMRG algorithm 
in the fitting interval (in our work T £ (0.05, 0.1)) is less 
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FIG. 5: The x cyc dependence of the spin gap for x = 0.5 and 
1.0 obtained by fitting the Eqs. Q (signed as C) and © 
(given as x), respectively. 
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FIG. 6: x/Xmax versus T/A for x = 0.5, 1.0 with a; cyc = 0, 
0.02, 0.04, 0.06, 0.08, and 0.1 in ascending order along the 
direction of arrows. 



than 10 -2 at most. On the other hand, in obtaining the 
fitting formulas ^ and ((7|), the low-temperature limit 
and quadratic approximation (See Eq. ([5])) have been 
adopted. For the specific heat, the expression including 
the higher order contributions of T/A is given explicitly 

aa31 



C(T) = 



3A 

1 - 
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T 

T 
A 



3/2 



-A/T 



(8) 
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The uncertainty in Eq. ((6]) is roughly estimated as 

2" 



SA ~ Tin 



1 



T 
A 



3 (T 

4 I A 



(9) 



where T represents some value of temperature T lying 
in the fitting interval. From this, in our case the uncer- 
tainty of gap due to the fitting formulas can be estimated 
approximately as SA « 0.01 — 0.02. Based on these con- 
siderations, we conclude that the error in the gap fitting 
is about 10 -2 , which is not exceeding the symbol size in 

Fig. m 





Mo 


Mi 


M2 


M3 


eo 


A 


Ac 





1.58 


0.59 


-0.38 


0.09 


-1.156 


0.53 


0.51 


0.02 


1.57 


0.60 


-0.39 


0.10 


-1.137 


0.47 


0.46 


0.04 


1.55 


0.62 


-0.41 


0.12 


-1.119 


0.42 


0.40 


0.06 


1.54 


0.63 


-0.42 


0.13 


-1.101 


0.37 


0.35 


0.08 


1.53 


0.64 


-0.43 


0.15 


-1.084 


0.31 


0.30 


0.1 


1.52 


0.65 


-0.44 


0.15 


-1.068 


0.27 


0.26 



TABLE I: The partition function fitting results by Eq. (fT0|) 
for the x = 1.0 case. Ac in the last column is the energy gap 
obtained by Eq. (|6). as plotted in Fig. [5] 



III. PARTITION FUNCTION AND THE 
DISPERSION RELATIONS 

In TMRG calculations, the partition function 
■Ztmrg(/3) is the original and the most accurate numer- 
ical quantity. In previous section we obtained the spe- 
cific heat and spin susceptibility, which can also be got 
from the derivatives of the partition function directly. In 
the following we will discuss the single-particle state in 
this system and show that its dispersion, as well as the 
ground-state energy and the magnitude of the gap, can 
be obtained from the analyses of the partition function. 

In the rung singlet phase, the elementary excitations 
are gapped. The energy gap and dispersion for the sin- 
gle particle state can be analyzed perturbatively^ When 
Jicg = Jcyc = 0, the system is consisted of an array of 
rungs with a vacuum state being the product of spin 
singlet state on each rung. Above this vacuum state, 
a hard core boson of spin 1 is the only possible excited 
state on each rung. The system has a flatband dispersion 
e(k) = J rung for single particle excitations. Due to the 
hard core scattering between the bosons, the statistics of 
the system is essentially Fermi-Dirac type^ 2 - which allows 
us to write down the partition function for the system: 

lnZ cS (P)= ln(l + 3e-P< k ^dk-/3e , (10) 

where eo is the ground-state energy per rung in the 
thermodynamic limit and the prefactor 3 before e~^ e ^ 
comes from the three internal states of the single particle 
excitations. 

In the rung singlet phase, switching upon the couplings 
Jieg and J cyc will induce a dispersion ^ 10 ' 20 for the single 
particle state, and the interaction between the particles. 
Classified by the internal symmetries, the two-particle 
interaction is cither spin-channel independent or spin- 
channel dependent. A bound state may also appear in 
certain spin-channels. This has been observed in a pre- 
vious study^ in which it is shown that the bound state 
is also momentum dependent. Therefore, the interaction 
is very complicated and depends on the total momentum 
of the scattering particles. 

Since we are dealing with a system with massive parti- 
cles at zero chemical potential, the density of the particles 



would be dilute at low temperatures. Taking into account 
of the form of the effective partition function for single 
particle dispersion Z e g, we expect it would match the ac- 
curate TMRG results at both high and low temperatures 
by including higher order harmonics in the dispersion^ 

e(fc) = mo + Mi cos(fc) + M2 cos(2fc) + M3 cos(3/c). (11) 

The energy gap is given by A = Mo — Ml + M2 — A*3 at 
k = 7r. Starting with a set of trial values for m's given 
by perturbative calculation^ we minimize the squared 
2-norm of the residual in free energy calculated from the 
trial dispersion: 

VF(Mo,Mi:M2,M3,e ) = 

N 

[ ln Z eff (ft) - hi Ztmrg (ft)] 2 , (12) 

i=l 

where ft are smoothly distributed in the interval (5, 15) 
with a spacing A/3 = 0.05. Typically, we show the fitting 
results for the x = 1.0 case in Table|TJ The corresponding 
dispersions are plotted in Fig. [7J In our final numerical 
results, the values of W in Eq. (fT2"|) are of order 10~ 6 , 
and the mean deviation at each data point is less than 
10" 4 . 

From Table [H we find that mo an d Mi monotonically 
decreases and increases with x cyc , respectively This is 
consistent with the perturbative results^ - On the other 
hand, the third order coefficient M3 increases with x cyc . 
At a; C y C = 0.1, it reaches 0.15, about one-third of M2- 
This indicates that the approximation up to the third 
order in the dispersion is necessary for x = 1.0 case. 
On the contrary, for x = 0.5, Mo ~ T Mi ~ 0.5 — 0.6, the 
absolute value of both M2 and M3 are less than 0.1, and the 
dispersion can be satisfactorily approximated by single- 
component harmonics with only Mi- This fact reflects 
that with the increase of J\ eg and J C y C , higher orders in 
the harmonic expansion are needed to describe the one- 
particle excitation in the spectrum. 

The dispersion relations shown in Fig. [7] should corre- 
spond to the one-triplet magnon branch in the spectrum 
of the ladder. We compare the results with the previous 
studies by the DMRGi 4 ^ We find they have quite similar 
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FIG. 7: (Color online) Dispersion (in units of Jr Ung ) for x — 
1.0 with x cyc = 0, 0.02, 0.04, 0.06, 0.08, and 0.1 in descending 
order along the direction of arrows. 



shapes and are quantitatively agreement in the vicinity 
of k = ir. The deviation when approaching k = can be 
attributed partially to the two-triplet bound state aris- 
ing from the magnon-magnon interaction that we have 
ignored^ 

The energy gaps obtained in this way are presented in 
Table HI and are compared with their counterparts from 
the specific heat data. We see that they are rather coin- 
cident. In addition, we have obtained the ground state 
energy (eo in Eq. (|10j) ). which is compatible with the 
results from the zero-temperature DMRG scheme. 



by the TMRG method. The temperature dependences of 
the specific heat and the susceptibility of the ladder in 
the rung-singlet phase are analyzed numerically. Based 
on the above results, we extracted the values of the spin 
gap and found that it manifests approximately linear be- 
havior with the cyclic interaction in the parameter space 
we considered. We also showed that the dispersion for 
low-lying excitations can be obtained from the partition 
function very effectively by combining with the mean- 
field or perturbative studies. Comparing with the high 
temperature series expansion and the exact diagonaliza- 
tion approaches we believe that the results presented 
in this study give more direct and accurate physical quan- 
tities in the thermodynamic limit. Further experimental 
results on specific heat and susceptibility may allow one 
to determine the various parameters for the two-leg spin 
ladder materials via this approach. In addition, the dis- 
persion relations extracted from the free energy may be 
used to be compared with the experimental data of the 
inelastic neutron scattering spectra. 
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